\chapter{Control Variable and State Transformations}
\label{ch:control_transforms}

The control variable framework in GSI provides the mathematical foundation for transforming between control variables (optimized in the minimization process) and physical state variables (used in the forward model). This chapter examines the implementation of these critical transformations and their role in variational data assimilation.

\section{Theoretical Foundation}

Control variables represent perturbations in a transformed space that facilitates optimization. The transformation system implements the relationship:

\begin{equation}
\mathbf{x} = \mathbf{x}_b + \mathbf{B}^{1/2} \hat{\mathbf{x}}
\label{eq:control_transform}
\end{equation}

where $\mathbf{x}$ is the state vector, $\mathbf{x}_b$ is the background state, $\mathbf{B}^{1/2}$ is the square root of the background error covariance matrix, and $\hat{\mathbf{x}}$ is the control vector.

The control variable parameterization includes:
\begin{itemize}
\item Stream function ($\psi$) and velocity potential ($\chi$) for horizontal wind
\item Unbalanced pressure and temperature components
\item Normalized relative humidity for moisture analysis
\item Surface pressure and other prognostic variables
\end{itemize}

\section{Control to State Transformation}

The \texttt{control2state} subroutine performs the forward transformation from control variables to physical state variables:

\subsection{Core Transformation Algorithm}

\begin{lstlisting}[language=Fortran,caption={Control to State Implementation},label=code:control2state]
subroutine control2state(xhat, sval, bval)
  use control_vectors, only: control_vector
  use gsi_bundlemod, only: gsi_bundle
  use bias_predictors, only: predictors
  
  type(control_vector), intent(in) :: xhat
  type(gsi_bundle), intent(inout) :: sval  
  type(predictors), intent(inout) :: bval
  
  ! Local variables
  real(r_kind), pointer, dimension(:,:,:) :: u_ptr, v_ptr
  real(r_kind), pointer, dimension(:,:,:) :: prs_ptr, tv_ptr, q_ptr
  real(r_kind), pointer, dimension(:,:,:) :: sf_ptr, vp_ptr
  real(r_kind), pointer, dimension(:,:) :: ps_ptr
  
  ! Check for sqrt(B) mode - not supported
  if (lsqrtb) then
    call die('control2state', 'not available in sqrt(B) mode')
  end if
  
  ! Variable discovery and pointer assignment
  call inquire_state_vars(sval, cw_exists, qi_exists, ql_exists)
  call gsi_bundlegetpointer(sval, 'u', u_ptr, istatus)
  call gsi_bundlegetpointer(sval, 'v', v_ptr, istatus)
  call gsi_bundlegetpointer(sval, 'tv', tv_ptr, istatus)
  call gsi_bundlegetpointer(sval, 'q', q_ptr, istatus)
  call gsi_bundlegetpointer(sval, 'ps', ps_ptr, istatus)
  
  ! Transform stream function and velocity potential to winds
  if (getindex(cvars3d, 'sf') > 0 .and. getindex(cvars3d, 'vp') > 0) then
    call gsi_bundlegetpointer(xhat, 'sf', sf_ptr, istatus)
    call gsi_bundlegetpointer(xhat, 'vp', vp_ptr, istatus)
    call getuv(sf_ptr, vp_ptr, u_ptr, v_ptr, 0)
  end if
  
  ! Transform pressure and temperature variables
  if (getindex(cvars3d, 'ps') > 0) then
    call gsi_bundlegetpointer(xhat, 'ps', ps_cv_ptr, istatus)
    ps_ptr = ps_cv_ptr
  end if
  
  ! Transform relative humidity to specific humidity
  if (getindex(cvars3d, 'q') > 0) then
    call gsi_bundlegetpointer(xhat, 'q', rh_ptr, istatus)
    call normal_rh_to_q(rh_ptr, tv_ptr, prs_ptr, q_ptr)
  end if
  
  ! Transform virtual temperature to sensible temperature
  if (getindex(svars3d, 'tsen') > 0) then
    call gsi_bundlegetpointer(sval, 'tsen', tsen_ptr, istatus)
    call tv_to_tsen(tv_ptr, q_ptr, tsen_ptr)
  end if
end subroutine control2state
\end{lstlisting}

\subsection{Wind Field Transformation}

The transformation from stream function and velocity potential to horizontal wind components uses the \texttt{getuv} subroutine:

\begin{equation}
\mathbf{u} = -\frac{\partial \psi}{\partial y} + \frac{\partial \chi}{\partial x}
\end{equation}
\begin{equation}
\mathbf{v} = \frac{\partial \psi}{\partial x} + \frac{\partial \chi}{\partial y}
\end{equation}

\begin{lstlisting}[language=Fortran,caption={Wind Transformation},label=code:wind_transform]
call getuv(sf_ptr, vp_ptr, u_ptr, v_ptr, 0)
! Computes horizontal derivatives using spectral transforms
! or finite difference methods depending on grid configuration
\end{lstlisting}

\subsection{Moisture Variable Transformation}

The moisture transformation converts normalized relative humidity back to specific humidity:

\begin{lstlisting}[language=Fortran,caption={Humidity Transformation},label=code:humidity_transform]
subroutine normal_rh_to_q(rh, tv, prs, q)
  ! Convert normalized RH to specific humidity
  ! rh: normalized relative humidity (control variable)  
  ! tv: virtual temperature
  ! prs: pressure
  ! q: specific humidity (state variable)
  
  do k = 1, nsig
    do j = 1, lon2  
      do i = 1, lat2
        ! Compute saturation mixing ratio
        call genqsat(qsat, tv(i,j,k), prs(i,j,k))
        
        ! Convert normalized RH to mixing ratio
        qmix = rh(i,j,k) * qsat / (1.0_r_kind + qsat)
        
        ! Convert mixing ratio to specific humidity
        q(i,j,k) = qmix / (1.0_r_kind + qmix)
        
        ! Apply bounds checking
        q(i,j,k) = max(qmin, min(qmax, q(i,j,k)))
      end do
    end do  
  end do
end subroutine normal_rh_to_q
\end{lstlisting}

\section{Adjoint Transformation}

The adjoint transformation \texttt{control2state\_ad} implements the transpose of the control-to-state operator:

\subsection{Adjoint Implementation}

\begin{lstlisting}[language=Fortran,caption={Adjoint Control to State},label=code:control2state_ad]
subroutine control2state_ad(xhat, sval, bval)
  ! Adjoint of control2state transformation
  ! Accumulates gradients from state space back to control space
  
  type(control_vector), intent(inout) :: xhat
  type(gsi_bundle), intent(in) :: sval
  type(predictors), intent(in) :: bval
  
  ! Adjoint of temperature transformation
  if (getindex(svars3d, 'tsen') > 0) then
    call gsi_bundlegetpointer(sval, 'tsen', tsen_ptr, istatus)
    call tv_to_tsen_ad(tv_ptr, q_ptr, tsen_ptr)
  end if
  
  ! Adjoint of humidity transformation
  if (getindex(cvars3d, 'q') > 0) then
    call gsi_bundlegetpointer(xhat, 'q', rh_ptr, istatus)  
    call normal_rh_to_q_ad(rh_ptr, tv_ptr, prs_ptr, q_ptr)
  end if
  
  ! Adjoint of wind transformation
  if (getindex(cvars3d, 'sf') > 0 .and. getindex(cvars3d, 'vp') > 0) then
    call gsi_bundlegetpointer(xhat, 'sf', sf_ptr, istatus)
    call gsi_bundlegetpointer(xhat, 'vp', vp_ptr, istatus)
    call getuv_ad(sf_ptr, vp_ptr, u_ptr, v_ptr, 0)
  end if
end subroutine control2state_ad
\end{lstlisting}

\subsection{Adjoint Code Verification}

The adjoint transformations satisfy the adjoint test condition:
\begin{equation}
\langle \mathbf{L} \mathbf{x}, \mathbf{y} \rangle = \langle \mathbf{x}, \mathbf{L}^T \mathbf{y} \rangle
\label{eq:adjoint_test}
\end{equation}

\begin{lstlisting}[language=Fortran,caption={Adjoint Test Implementation},label=code:adjoint_test]
subroutine test_control2state_adjoint()
  ! Test adjoint accuracy
  real(r_kind) :: dot1, dot2, relative_error
  
  ! Initialize random control vector
  call random_cv(xhat1)
  call random_state(sval2)
  
  ! Forward transformation
  call control2state(xhat1, sval1, bval1)
  dot1 = dot_product(sval1, sval2)
  
  ! Adjoint transformation  
  call control2state_ad(xhat2, sval2, bval2)
  dot2 = dot_product(xhat1, xhat2)
  
  ! Verify adjoint accuracy
  relative_error = abs(dot1 - dot2) / max(abs(dot1), abs(dot2))
  
  if (relative_error > adjoint_tolerance) then
    write(*,*) 'Adjoint test failed:', relative_error
  else
    write(*,*) 'Adjoint test passed:', relative_error  
  end if
end subroutine test_control2state_adjoint
\end{lstlisting}

\section{Control to Model Transformation}

The \texttt{control2model} subroutine transforms control variables to model grid variables, handling interpolation between different resolutions:

\subsection{Multi-Resolution Support}

\begin{lstlisting}[language=Fortran,caption={Control to Model Transformation},label=code:control2model]
subroutine control2model(xhat, wbundle, bval)
  ! Transform control variables to model grid
  ! Supports dual-resolution configurations
  
  type(control_vector), intent(in) :: xhat
  type(gsi_bundle), intent(inout) :: wbundle
  type(predictors), intent(inout) :: bval
  
  ! Work arrays for grid transformations
  real(r_kind), allocatable :: work_cv(:,:,:)
  real(r_kind), allocatable :: work_model(:,:,:)
  
  ! First transform to state variables
  call control2state(xhat, sval, bval)
  
  ! Handle dual resolution if configured
  if (dual_res) then
    ! Interpolate from analysis grid to model grid
    call anl2model_interp(sval, wbundle)
  else
    ! Direct copy for same resolution
    wbundle = sval
  end if
  
  ! Apply balance constraints if configured
  if (l_balance) then
    call apply_balance_constraint(wbundle)
  end if
end subroutine control2model
\end{lstlisting}

\subsection{Balance Constraint Application}

Balance relationships maintain dynamical consistency:

\begin{lstlisting}[language=Fortran,caption={Balance Constraint Implementation},label=code:balance_constraint]
subroutine apply_balance_constraint(wbundle)
  ! Apply geostrophic and hydrostatic balance
  
  ! Geostrophic balance: f * u = -∂Φ/∂y, f * v = ∂Φ/∂x
  if (l_geostrophic_balance) then
    call compute_geostrophic_balance(u_ptr, v_ptr, phi_ptr)
  end if
  
  ! Hydrostatic balance: ∂Φ/∂p = -RT/p
  if (l_hydrostatic_balance) then
    call compute_hydrostatic_balance(phi_ptr, t_ptr, p_ptr)
  end if
  
  ! Nonlinear normal mode constraint
  if (l_strong_constraint) then
    call strong_baldiag_inc(wbundle)
  end if
end subroutine apply_balance_constraint
\end{lstlisting}

\section{Ensemble Control Variable Framework}

The hybrid ensemble-variational system requires specialized transformations for ensemble control variables:

\subsection{Ensemble Control to State}

\begin{lstlisting}[language=Fortran,caption={Ensemble Control Transformation},label=code:ensctl2state]
subroutine ensctl2state(xhat_ens, sval, bval)
  ! Transform ensemble control variables to state
  
  type(control_vector), intent(in) :: xhat_ens
  type(gsi_bundle), intent(inout) :: sval
  type(predictors), intent(inout) :: bval
  
  ! Ensemble perturbation reconstruction
  if (l_hyb_ens) then
    ! Standard ensemble transform
    call ensemble_forward_model(xhat_ens, sval_ens)
    
    ! Apply localization if configured
    if (l_loc) then
      call apply_localization(sval_ens, sval)
    else
      sval = sval_ens
    end if
    
  end if
  
  ! Handle dual resolution ensemble
  if (dual_res .and. l_hyb_ens) then
    call ens2anl_ad_bilin(sval_ens_low, sval_high)
  end if
end subroutine ensctl2state
\end{lstlisting}

\subsection{Localization Implementation}

Ensemble localization limits spurious correlations:

\begin{lstlisting}[language=Fortran,caption={Ensemble Localization},label=code:localization]
subroutine apply_localization(sval_ens, sval_loc)
  ! Apply Gaspari-Cohn localization function
  
  real(r_kind) :: distance, correlation
  real(r_kind), parameter :: localization_radius = 500.0e3_r_kind
  
  do k = 1, nsig
    do j = 1, lon2
      do i = 1, lat2
        ! Compute distance from localization center
        distance = great_circle_distance(lat(i), lon(j), &
                                       lat_center, lon_center)
        
        ! Apply Gaspari-Cohn function
        if (distance <= localization_radius) then
          correlation = gaspari_cohn_function(distance/localization_radius)
          sval_loc(i,j,k) = correlation * sval_ens(i,j,k)
        else
          sval_loc(i,j,k) = 0.0_r_kind
        end if
      end do
    end do
  end do
end subroutine apply_localization

real(r_kind) function gaspari_cohn_function(r)
  ! Fifth-order piecewise rational function
  real(r_kind), intent(in) :: r
  
  if (r <= 1.0_r_kind) then
    gaspari_cohn_function = 1.0_r_kind - (5.0_r_kind/3.0_r_kind)*r**2 + &
                           (5.0_r_kind/8.0_r_kind)*r**3 + &
                           (1.0_r_kind/2.0_r_kind)*r**4 - &
                           (1.0_r_kind/4.0_r_kind)*r**5
  else if (r <= 2.0_r_kind) then
    gaspari_cohn_function = 4.0_r_kind - 5.0_r_kind*r + &
                           (5.0_r_kind/3.0_r_kind)*r**2 + &
                           (5.0_r_kind/8.0_r_kind)*r**3 - &
                           (1.0_r_kind/2.0_r_kind)*r**4 + &
                           (1.0_r_kind/12.0_r_kind)*r**5 - &
                           2.0_r_kind/(3.0_r_kind*r)
  else
    gaspari_cohn_function = 0.0_r_kind
  end if
end function gaspari_cohn_function
\end{lstlisting}

\section{Cloud and Hydrometeor Transformations}

The GSI system includes specialized transformations for cloud and hydrometeor variables:

\subsection{Cloud Water Transformation}

\begin{lstlisting}[language=Fortran,caption={Cloud Water Transformation},label=code:cw_transform]
subroutine cw2hydro_tl(cw, ql, qi, qr, qs, qg, qh)
  ! Transform cloud water mixing ratio to hydrometeor species
  ! Uses Thompson or other microphysics scheme
  
  real(r_kind), intent(in) :: cw(:,:,:)
  real(r_kind), intent(out) :: ql(:,:,:), qi(:,:,:)
  real(r_kind), intent(out) :: qr(:,:,:), qs(:,:,:)
  real(r_kind), intent(out) :: qg(:,:,:), qh(:,:,:)
  
  real(r_kind) :: temperature, pressure
  real(r_kind) :: cw_total, cw_liquid, cw_ice
  
  do k = 1, nsig
    do j = 1, lon2
      do i = 1, lat2
        temperature = ges_tv(i,j,k)
        pressure = ges_prsl(i,j,k)
        cw_total = cw(i,j,k)
        
        ! Partition based on temperature
        if (temperature > t_liquid) then
          ! All liquid above freezing
          cw_liquid = cw_total
          cw_ice = 0.0_r_kind
        else if (temperature < t_ice) then
          ! All ice below freezing  
          cw_liquid = 0.0_r_kind
          cw_ice = cw_total
        else
          ! Mixed phase - linear interpolation
          alpha = (temperature - t_ice) / (t_liquid - t_ice)
          cw_liquid = alpha * cw_total
          cw_ice = (1.0_r_kind - alpha) * cw_total
        end if
        
        ! Assign to output arrays
        ql(i,j,k) = cw_liquid
        qi(i,j,k) = cw_ice
        
        ! Set precipitation species to zero (not directly controlled)
        qr(i,j,k) = 0.0_r_kind
        qs(i,j,k) = 0.0_r_kind
        qg(i,j,k) = 0.0_r_kind
        qh(i,j,k) = 0.0_r_kind
      end do
    end do
  end do
end subroutine cw2hydro_tl
\end{lstlisting}

\section{Memory Management and Performance}

The transformation system implements efficient memory management:

\subsection{Temporary Array Management}

\begin{lstlisting}[language=Fortran,caption={Memory Efficient Transformations},label=code:memory_efficient]
subroutine efficient_control2state(xhat, sval, bval)
  ! Use temporary work arrays to minimize memory footprint
  
  type(gsi_bundle) :: work_bundle
  real(r_kind), allocatable :: work_3d(:,:,:)
  
  ! Allocate minimal work space
  call gsi_bundlecreate(work_bundle, grid, 'work_bundle')
  allocate(work_3d(lat2,lon2,nsig))
  
  ! Process variables sequentially to minimize peak memory
  call process_wind_variables(xhat, sval, work_bundle)
  call process_mass_variables(xhat, sval, work_bundle)  
  call process_moisture_variables(xhat, sval, work_bundle)
  call process_surface_variables(xhat, sval, work_bundle)
  
  ! Clean up
  deallocate(work_3d)
  call gsi_bundledestroy(work_bundle)
end subroutine efficient_control2state
\end{lstlisting}

\subsection{Parallel Processing Optimization}

\begin{lstlisting}[language=Fortran,caption={Parallel Transformation},label=code:parallel_transform]
! OpenMP parallel sections for independent transformations
!$OMP PARALLEL SECTIONS
!$OMP SECTION
  ! Wind transformation
  call getuv(sf_ptr, vp_ptr, u_ptr, v_ptr, 0)
!$OMP SECTION  
  ! Mass variable transformations
  call getprs_tl(ps_ptr, tv_ptr, prsl_ptr)
!$OMP SECTION
  ! Moisture transformations  
  call normal_rh_to_q(rh_ptr, tv_ptr, prsl_ptr, q_ptr)
!$OMP END PARALLEL SECTIONS
\end{lstlisting}

\section{Quality Control Integration}

The transformation system integrates with quality control mechanisms:

\subsection{Physical Constraint Application}

\begin{lstlisting}[language=Fortran,caption={Physical Constraints},label=code:physical_constraints]
subroutine apply_physical_constraints(sval)
  ! Apply physical bounds and constraints
  
  ! Humidity constraints
  do k = 1, nsig
    do j = 1, lon2
      do i = 1, lat2
        ! Prevent supersaturation
        if (clip_supersaturation) then
          call genqsat(qsat, tv_ptr(i,j,k), prsl_ptr(i,j,k))
          q_ptr(i,j,k) = min(q_ptr(i,j,k), factqmax * qsat)
        end if
        
        ! Minimum humidity
        q_ptr(i,j,k) = max(q_ptr(i,j,k), qmin)
        
        ! Temperature bounds
        tv_ptr(i,j,k) = max(tmin, min(tmax, tv_ptr(i,j,k)))
      end do
    end do
  end do
  
  ! Surface pressure bounds
  do j = 1, lon2  
    do i = 1, lat2
      ps_ptr(i,j) = max(psmin, min(psmax, ps_ptr(i,j)))
    end do
  end do
end subroutine apply_physical_constraints
\end{lstlisting}

\section{Summary}

The control variable and state transformation framework in GSI provides:

\begin{itemize}
\item \textbf{Mathematical rigor}: Proper implementation of control variable theory
\item \textbf{Physical consistency}: Maintenance of dynamical balance relationships  
\item \textbf{Ensemble integration}: Seamless hybrid ensemble-variational capability
\item \textbf{Multi-resolution support}: Efficient dual-resolution configurations
\item \textbf{Adjoint accuracy}: Precise adjoint implementations for gradient computation
\item \textbf{Performance optimization}: Memory-efficient and parallel-optimized code
\item \textbf{Quality control}: Integrated physical constraint enforcement
\end{itemize}

These transformations form the core of the GSI variational data assimilation system, enabling optimal blending of observations with background forecasts while maintaining physical consistency and computational efficiency.